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Abstract 



A scheme is presented to decompose the exponential of a two-body opera- 
tor in a discrete sum over exponentials of one-body operators. This discrete 
decomposition can be used instead of the Hubbard-Stratonovich transforma- 
tion in auxiliary-field quantum Monte-Carlo methods. As an illustration, the 
decomposition is applied to the Hubbard model, where it is equivalent to the 
discrete Hubbard-Stratonovich transformation introduced by Hirsch, and to 
the nuclear pairing Hamiltonian. 
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I. INTRODUCTION 



In auxiliary- field quantum Monte-Carlo methods (AFQMC), such as the projector, 
grand-canonical |IJ and shell-model quantum Monte-Carlo methods @, the Boltzmann op- 
erator e _/3h , with h the Hamiltonian, is decomposed in a sum or integral of exponentials of 
one-body operators. This sum or integral is then evaluated using Monte-Carlo techniques. 
For the decomposition, these methods rely on the Hubbard-Stratonovich transformation 
0.0], which is based on the identity 

e /3p 2 = _L r e-^e^da, (1) 

V 2,TT J —oo 

where p is a one-body operator. In order to avoid problems due to non-commuting oper- 
ators, one can split up the Boltzmann operator using the Suzuki- Trotter formula [||]. One 
can discretize the Hubbard-Stratonovich transformation by applying a Gaussian quadrature 
formula to the integral over a. After a Suzuki- Trotter expansion in N t slices, a three-points 
quadrature formula leads to an error of the order of O (|U 3 ). This is of the same order 
in (3 and N t as the errors due to the non-commutativity of the squared operators that build 
up h. For some systems one can derive an exact, discrete Hubbard-Stratonovich transform. 
Hirsch showed how one can write an operator of the form e~ l3Un ^ ni , where n a is the site 
occupation number for an electron with spin projection a, exactly as a sum of two exponen- 
tials of one-body operators 0. Recently, Gunnarsson and Koch extended this to systems 
with higher orbital degeneracy ]/]]. 

The aim of this paper is to describe another discrete decomposition scheme, which is exact 
for a certain class of operators. This decomposition scheme is generalized to any two-body 
Hamiltonian using the Suzuki- Trotter formula. For the application in AFQMC methods, 
especially the shell-model Quantum Monte-Carlo method, this new decomposition has the 
advantage, compared to the discretized Hubbard-Stratonovich transform based on Eq. ([[]), 
that it is more accurate and that it leads to low-rank matrices. This leads to faster matrix 
multiplications and requires less computer memory. AFQMC methods for fermions often 
have sign problems [IJ . Fahy and Hamann || showed that these sign problems can be related 
to the diffusive behavior of states in the Hubbard-Stratonovitch transformation. Because our 
decomposition, in general, is based on exponentials of one-body operators of a completely 
different type, one can expect different sign properties. Our decomposition is not free of 
sign problems, but there might be systems where it leads to a sign rule while the Hubbard- 
Stratonovich transformation does not, or where our decomposition causes significally less 
sign problems. 

In Section II we introduce a matrix notation for Slater determinants and operators needed 
for a clear discussion of the decomposition. In Section III a basic lemma is given on which 
the decomposition is based. In Section IV the exact decomposition for a certain class of 
operators is presented. We indicate how to apply this decomposition to a general two-body 
Hamiltonian. In Section V the relation with Hirsch's decomposition for the Hubbard model 
is elucidated. Finally, in Section VI the decomposition for the nuclear pairing Hamiltonian 
is discussed and illustrated with AFQMC-results for an exactly solvable model. 
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II. A MATRIX NOTATION FOR SLATER DETERMINANTS AND OPERATORS 



In order to avoid confusion between matrix representations in the space of single-particle 
states and the operators themselves in Fock space, we will denote the former with upper case 
and the latter with lower case characters. Let {<pi, . . . , 4>n s } be a set of basis states for the 
single-particle space, a±, . . . ,cin s be the related creation operators and the A-particle state 
the antisymmetrized product of a set of single-particle states Mijfa, j = 1 ... A. 
i.e. ^ m is a Slater determinant. Thus in second quantization one can write 

A /N s \ 

i*> = n e^a i>- ( 2 ) 

j= i \i=i / 

This defines a matrix representation M for a Slater determinant ^m- The value of this 
representation is that one can represent certain operations on the Slater determinant by 
matrix operations on M. e.g. the overlap between two Slater determinants ^ m 1 and ^ m is 
given by (^m\^M') — det (^M^M'^j. The exponential of a one-body operator acting on 

results in a new Slater determinant, e'^^u) = \^?M') (this is a corollary of the Thouless- 
theorem U), whose matrix representation is related to M by M' = e~ l3H M, where the 
N s x N B matrix H is defined by h = HijCqaj. Reversily, given a N s x N s matix Q, one 
can consider the operator 0(Q), defined by its action on Slater determinants: 

6{Q) : |^ M ) — > 6{Q)\* M ) = |*Af/) with M' = QM. (3) 

If Q is non-singular, 0(Q) is the exponential of a one-body operator. 



III. A BASIC LEMMA 

Lemma: The operation represented by the unit matrix plus a matrix of rank two can 
be expressed as a sum of one- and two-body operators in the following way: 

6{l + a B\B 4 + (3 B\B Z ) = 1 + a b% + (3 S|S 3 + a(3 b\ti 2 hk, (4) 

where B\,B 2 ,B 3 and B± are 1 x N s row matrices and b\ = Y,j=i{Bk)j , k = 1,2,3,4. 
Proof: Consider the A-particle Slater determinant ^ m represented by the matrix M. 
Consider also a Slater determinant ^l, that has particles in the single-particle states 
4>i 2 i ■ ■ ■ ) a- The Slater determinants of this type constitue a basis of the A-particle 
Hilbert space. The overlap of ^ l with ^ m is given by 

(* £ |*iif> = det (M.x M.2 ■ ■ • M, A ) . (5) 

The notation M_j denotes the vector that is given by the j th column of M, the notation 
B for an iV-element vector B denotes the A-element vector (B^B^ ■ ■ ■ Bi A ). The operator 
in Eq.(|) transforms * M into * M / with M' = (1 + aB{B 4 + f3B ] 2 B 3 )M. To calculate the 
overlap of *&m' with ty L , we have to replace every column Mj in Eq.(H): 

M_j — ► M'j — M_j + atjB\ + PjBl, with a s = a B 4 M,j, (3j — (3 B 3 M_j. (6) 
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The overlap is then given by 

(* L \* M >) = det (M A + a x B\ + (5 X B\ ■■■ M. A + a A B{ + /3 A Bl) . (7) 

This determinant can be expanded as the sum of all determinants that are obtained by 
selecting in every column of Eq.(^) one of the terms M_j, ajB\ or l3jB\. If in more than one 
column the term 0ijB x is selected, then the determinant has two linearly dependent columns, 
so it will vanish. The same holds for the term (3jB\. Only four types of determinants remain: 

• M is selected in every column. This determinant is just {^l\^m) ( see Eq.(f|)). 

• ajB{ is selected in column j, M in all others. These determinants sum up to 

L \ab\i>4\^> M ) (one particle is moved from state 6 4 to state b x ). 

• (3jB\ is selected in column j, M in all others. These determinants sum up to 

LlPb^bsl^ m) (one particle is moved from state 63 to state b-z). 

• ctjB\ is selected in column j, fikBl is selected in column k, M, in all others. These 
determinants sum up to (\& ila/^Jfol^M 1 ^ m) (two particles are moved from states 64 
and 63 to states b\ and b?). 

Taking all these terms together, we find that 

(V l \Vm>) = 1 + ab\b 4 + pQb 3 + a/3b\blhh\^ M ) ■ (8) 
This holds for any basis state so that 

*Af = (l + ab\h + (3blh + apPMsh) * M . (9) 

This proves Eq.(£|). 

End of proof. 

IV. A DISCRETE HUBABRD STRATONOVICH DECOMPOSITION 

Consider a fermionic two-body operator Q of the form 

q = Qij ( B i)k ( B 2)i ala]aia k . (10) 

i,j,k,l=l 

An operator of this form has the special property that 

q 2 = Xq, with X= J2 (Qki ~ Qik) {B l ) k (B 2 ) l . (11) 
k,l=l 

Because of this relation, the exponential of q can be written as 
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Now we can use the lemma to obtain a discrete decomposition of e " q in a sum of exponentials 
of one-body operators: 

N " 1 I O ■ • I 

e -/3<? = ^ - ^ — (l + axijalh + ax' i:j a}b2 + XijX^a\a)b 2 b^ , 

ij=l A er=-l,+l W 

E E W^^ + ^A^i + ^t^), (13) 
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with the 1 x iV s row matrix which has a 1 on the i entry and zeros anywhere else, and 

e=EI3d, (i4) 

x<* = v^®> ( 15 ) 

X<< = si g n (7%) > (16) 



6 fc = £(£ fe )^, fc = l,2. (17) 

1=1 

This is an exact Hubbard- Stratonovich-like decomposition of the form of Eq. fli~0|) , To apply 
this discrete Hubbard-Stratonovich-like decomposition to the Boltzmann operator e~^ v for 
a general fermionic two-body operator v, one has to rewrite v as a sum of operators of the 
form Eq. (|T0|) . A trivial way to do this, is given by 

N s N s I N s \ 

v = V ijkl a\a]aia k = Qta, with q u = ^ V ijki a\a] a ; a fc . (18) 

i,j,k,l=l k,l=l \*J=1 / 

The Suzuki- Trotter formula can be used to split up the Boltzmann operator into factors 
with only one operator q k \ in the exponent. Then each of these factors can be decomposed 
exactly using the discrete decomposition in Eq. (|i3|) . Note that the total decomposition is 
no longer exact because of the non-commutativity of the operators q k \. The error will be of 
the order (D(P 3 /N t 2 ). It will be much smaller than in case of a decomposition based on a 
Gaussian discretization of the integral in Eq.(|]), because now the error is proportional to 
the commutators of the operators q k i and no longer to a power of v. 



V. RELATION TO HIRSCH'S DECOMPOSITION FOR THE HUBBARD 

HAMILTONIAN 

For the Hubbard model we have to find a decomposition for a Boltzmann operator of 
the form e~ /3Un ^ n - 1 . where U is the interaction strength and = a] J a ( j. a =f, [ is an index 
for the spin degree-of-freedom. The exponent has a two-body operator n^n^, which is an 
operator of the form of Eq.flToD, so we can apply the decomposition given in Eq. flilf ) and 
obtain: 

e -WH = l £ 6 (1 + axM + axM) , (19) 

<T=-1,+1 
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with N^(Ni) the matrix which is zero everywhere except for the diagonal element related to 
the spin-up (spin-down) site, which is equal to 1. x an d x' are given by 

X T = -Xi = VI ~ e-W for (3U > 0, (20) 
or XT = Xl = Ve-? u - 1 for pU < 0. (21) 

Now one could scale each term in Eq. fli~9"l) with an operator of the form e _/3/i ( nT+ni ) ; because 
in the canonical ensemble this is just a constant. The matrices in the decomposition now 
have to be multiplied with the matrix 1 + (e~^ — l) iV-j- + {eT^^ — \\ JVj_. In case of the 
repulsive Hubbard model, the choice \i = —U/2 leads to the discrete Hubbard- Stratonovich 
transform of Hirsch ||. From the computational point of view this particular choice of \i 
has the advantage that the matrix representation for the spin-down part is related to the 
matrix representation for the spin-up part by a matrix inversion. Then one only has to keep 
track of the spin-up part in actual AFQMC calculations. Hirsch's decomposition for the 
attractive Hubbard model can also be obtained from Eq.flT9D, with a particular choice for /i. 
In this case however, there is no computational advantage in taking any particular value. 



VI. APPLICATION TO THE NUCLEAR PAIRING HAMILTONIAN 

The Hamiltonian for nuclear pairing in a degenerate shell is given by 

N S 



k,k'=l 



h = -G a\a\a- k ,a kl . (22) 



Here it is assumed that there are 2Ng degenerate single-particle states. The single-particle 
energy is shifted to MeV. So there is no one-body part in the Hamiltonian. The states 
with j z > are labeled from 1 to Ng and k denotes the time-reversed state of state k. The 
many-body problem for this model can be solved analytically using the seniority scheme 

Using the Suzuki- Trotter formula, the Boltmann operator for this Hamiltonian can be 
written as 

e -ph = e -f<h e -f h . . . e -§4N S e -^N S . . . e -f& e -§& + q ^ ; (23) 

with 

I n s \ 

Qk = ~G ( J] a\,a\, \ a- k a k . (24) 

The error is of the order O (/3 3 ). It is assumed that (3 is small. In practice, one has to split 
(3 in a number of inverse-temperature slices using the Suzuki- Trotter formula. Then one can 
apply the procedure that is discussed here to each inverse-temperature slice seperately. We 
have ql = —Gq k . So we can find a decomposotion of the type given in Eq.(pT3|) 



2, 



e 



1 

E E w -6(l + a X AlA k + axAlA i ), (25) 



fc'=l CT =-l, +i 2N s 
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where \ 2 — As (e^ 2 ' ~ l) • 

We have applied this decomposition to study a degenerate shell of 20 states (N s = 10), 
with 10 particles. This could model the valence model space for neutrons in the fp shell 
in 56 Fe, if one neglects the energy gap between the lfr and the 2p3, If 5, 2pi orbitals. 
For the strength of the interaction we took G = 20 /A MeV = 20/56 MeV, as recom- 
mended in |jTI}| . We have performed a shell-model quantum Monte-Carlo calculation in 
the canonical ensemble, following 0, but now using the new decomposition of Eq . (p5[) in- 
stead of the Hubbard- Stratonovich transformation. In order to make the systematic error 
smaller than the statistical error, the inverse temperature (3 was split into slices of length 
0.05 MeV" 1 . We point out that the form 1 + a\A\,A k + axA\,A k can be rewritten as 

(l + axA\,A])j (l + axA^Aj^j, such that there is a symmetry between states with j z > 
and their time-reversed states. This symmetry guarantees that there will be no sign problem 
for systems with an even number of particles. This sign-rule is analogous to the sign rule 
for the pairing-plus-quadrupole Hamiltonian decomposed using the Hubbard- Stratonovich 
transform | I2| . In figure [l| we show the internal energy of the system as a function of tem- 



perature. In figure |2] we show the corresponding specific heat of the system as a function 
of temperature. The Monte-Carlo results are in excellent agreement with the analytical 
results. The peak in de specific-heat curve around a temperature of 0.8 MeV. can be as- 
sociated with the breakup of .P = + pairs. It is straightforward to take into account the 
different single-particles energies and more general forms of the pairing Hamiltonian: 

N S N s 

h = ^e k (h k + h k ) - G k ^a\a\aya k '- (26) 

fc=l k,k'=l 

Extension to even more general two-body Hamiltonians is possible. Then there can be sign 
problems at low temperatures. 



VII. CONCLUSION 

We have presented a new type of discrete Hubbard- Stratonovich decomposition for the 
Boltzmann operator. It is exact for a special class of two-body operators. Applied to the 
Hubbard Hamiltonian, it leads to Hisrch's discrete Hubbard-Stratonovich decomposition. 
The decomposition is well suited for the nuclear pairing Hamiltonian, where it leads to a 
sign rule for systems with an even number of particles. Quantum Monte-Carlo results based 
on this decomposition are in excellent agreement with the analytical results for an exactly 
solvable model. 
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FIGURES 
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FIG. 1. Internal energy versus temperature. Error bars on the Monte-Carlo data are omitted 
because they are smaller than the symbols marking the data points. 
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FIG. 2. Specific heat versus temperature. Error bars on the Monte-Carlo data represent 
95%-confidence intervals. 
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